#ifndef SCHRODINGEREQUATION_H_INCLUDED
#define SCHRODINGEREQUATION_H_INCLUDED

#include <utility>
#include <cmath>
#include <complex>
#include <functional>

#include "TypePoint.h"

namespace SchrodingerEquation{

constexpr double a0 = 0.02;

inline uint64_t GetFactorial(int k)
{
    static uint64_t memo[15];
    if(k >= 15)
        throw std::domain_error("Trying to calculate a factorial that will cause overflow.");
    if(k <= 1)
        return 1ULL;
    else
        return (memo[k] == 0 ? memo[k] = GetFactorial(k-1)*k : memo[k]);
}

inline int GetC(int m, int n)
{
    static uint64_t memo[15][15];

    if(m >= 15 || n >= 15)
        throw std::domain_error("Trying to calculate a C too large.");

    if(memo[m][n] == 0)
    {
        memo[m][n] = GetFactorial(n) / GetFactorial(m) / GetFactorial(n - m);
    }

    return memo[m][n];
}

inline double GetLaguerre(int a, int n, double x)
{
    double result = 0;
    for(int k = 0; k <= n; k++)
        result += std::pow(-1, k) * GetC(n - k, n + a) * std::pow(x, k) / GetFactorial(k);
    return result;
}

inline double GetLegendre(int n, double x)
{
    if(n == 0)
        return 1;
    if(n == 1)
        return x;
    return ((2.0 * n - 1) / n) * GetLegendre(n-1,x) - (n / (n+1)) * GetLegendre(n-2, x);
}

inline double GetAssocLegendre(int m, int l, double x)
{
    if(m == 0)
        return GetLegendre(l, x);
    return (x * (2 * l - 1) * GetAssocLegendre(m, l-1, x) - (l + m - 1) * GetAssocLegendre(m, l-2, x)) / (l-m);
}

inline std::complex <double> GetExpImag(double coef)
{
    // e^ix = cos(x) + i * sin(x)
    return std::complex<double>(std::cos(coef), std::sin(coef));
}

inline std::complex<double> HydrogenGenericComplex (std::pair<double, double> p, int n, int l, int m = 0)
{
    double Z = 1;
    const double & r = p.first;
    const double & theta = p.second;
    const double & phi = 0;
    double coeffR = 0;
    double coeffY = 0;

    coeffR = -std::sqrt(std::pow(2 * Z / n / a0, 3) * GetFactorial(n - l - 1) / (2 * n * GetFactorial(n + l)));
    coeffY = std::sqrt((2.0 * l + 1) / (4 * M_PI) * GetFactorial(l - std::abs(m)) / GetFactorial(l + std::abs(m)));

    std::complex<double> R, Y;

    R.real(coeffR * std::exp(-Z * r / n / a0) * std::pow(2 * Z * r / n / a0, l) * \
           GetLaguerre(2*l+1, n-l-1, 2 * Z * r / n / a0));
    Y.real(coeffY * GetAssocLegendre(m, l, std::cos(theta)));
    Y *= GetExpImag(m * phi);
    return (-R) * Y;
}

inline double HydrogenGeneric(std::pair<double, double> p, int n, int l, int m = 0)
{
    return std::norm(HydrogenGenericComplex(p, n, l, m));
}

};

using SchrodingerEquation::HydrogenGeneric;

#endif // SCHRODINGEREQUATION_H_INCLUDED
